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Consider the Dirichlet-to-Neumann operator Af in the exterior problem for 
the 2D Helmholtz equation outside a bounded domain with smooth boundary. 
Using parametrization of the boundary by normalized arclength, we treat Af as a 
pseudodifferential operator on the unit circle. We study its discrete symbol. 

We put forward a conjecture on the universal behaviour, independent of shape 
and curvature of the boundary of the symbol as the wavenumber k — > oo. The 
conjecture is motivated by an explicit formula for circular boundary and confirmed 
numerically for other shapes. It also agrees, on a physical level of rigor, with 
KirchhofFs approximation. The conjecture, if true, opens new ways in numerical 
analysis of diffraction in the range of moderately high frequencies. 

Introduction 

This work is a part of research aimed at an accurate and robust numerical algorithm for 
diffraction problems in mid-high frequency range, where the standard Boundary Inte- 
gral Equation methods fail due to large matrix size and, more importantly, to numerical 
contamination in quadratures. A natural idea to use the knowledge of geometric phase 
and to separate fast oscillations from slowly varying amplitudes has been converted to 
a practical method EQ,[2] with recent enhancements [7j. A drawback of that approach 
occurs in the presence of flattening boundary regions, where Kirchhoff 's amplitude be- 
comes singular. From numerical analyst's point of view, a method that has problem 
with small curvature is anti-intuitive. 

The point of our approach is to look for an object in theory whose high- 
frequency asymptotics stands flattening well and isn't sensitive to convexity 
assumptions. We suggest that the symbol of Dirichlet-to-Neumann operator might 
be such an object. 

We consider the 2D case and don't claim a ready-made extension of our results in 
3D. As a technical reason, we need a well defined full symbol of a pseudodifferential 
operator on a compact manifold (the boundary). In the 2D case, the boundary is a 
closed curve, so a special version of the PDO theory with discrete frequency variable is 
applicable, which attends to smooth kernels and doesn't require partitions of unity. 



Dirichlet-to-Neumann operator 



Let f2 be a bounded domain in R 2 with smooth boundary T. The exterior Dirichlet 
problem for the Helmholtz equation in polar coordinates r, <f> reads 
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d 2 u 1 du 1 d 2 u , 9 9 - 

Au = in + -ir + ^in2 = ~ k u m R \ n > 

^ M -7 / -l/2\ (1) 

— — — iku = o(r 'J, r — > oo, 
ar 

«|r = /• 

For a function / 6 C 1 (r), the problem has unique solution u, and its normal derivative 
g = <9 ra u|r is a continuous function on V. (For sharper conditions see e.g. [JJ], ^0] •) The 
map AA : / — ► g is called the Dirichlet-to- Neumann operator for Problem (|Tj). 

Operator A/" for the exterior of the unit disc 

Here the boundary is the unit circle S and it can be parametrized by (p. Consider 
Fourier series of the 2-7r-periodic functions /(</>) = u\s and g{4>) = d n u\s 

m = E KnY n \ g(<t>) = E 9{n)e m t 

n£Z ngZ 

The Helmholtz equation has outgoing elementary solutions in the product form 

where H$ are Hankel functions [B] (9.1.3). The solution u{r,4>) can be represented as 
a linear combination of the elementary solutions. Matching Fourier coefficients in the 
boundary data, we find 

g(n) = cr(n; k)f(n), 

where (cf. [S] (9.1.27.4) ) 



a(n; k) = k = -k " ^ + \n\. (2) 

The operator AA can be written in a pseudo differential fashion, the function a being 
its discrete symbol (here dependence of a on the wavenumber k is irrelevant and is 
omitted) 

AW) = E ff(n)/(n)e in * (3) 

neZ 

ASYMPTOTICS OF THE SYMBOL 

For a fixed and n — ► 00, we derive from [B] (9.3.1) 

#P(fc) ' 
so (in full agreement with pseudodifferential calculus) 

cr(n; fe) ~ — \n\, \n\ —> 00. (4) 
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On the other hand, if n is fixed, then by 8 u (9.2.3) 

<r(n; k) ~ ik, k — > oo. (5) 

The next result, in which the ratio t = n/k is fixed, interpolates between the above two 
special cases. Since n is an integer, the integral part function [•] is involved. 

Lemma. For any fixed t > and n = n(k,t) = [kt], 



n; fe) . , A J *Vl-t 2 , if t < 1 
hm — ^ — = oito t = { (6) 

This fact can be derived laboriously using ^8 (9.3.37-46) and asymptotics for the Airy 
functions. Instead, we demonstrate a simple argument, which quickly produces the 
formula in the case t ^ 1, and can be converted to a formal proof. 

Consider a recurrence for the ratios \i v of Hankel functions of orders v + 1 and v 
with fixed argument k. According to |8] (9.1.27.1), we have 

Mv + AV-i = 2v/k, 

or in the explicit difference form, 

fhs - fMs-1 = -fJv-1 - HZ-l + 2 ^> ^ = v / k - ( 7 ) 

The ratio t v varies slowly. Consider the difference equation with v ~ n and frozen 
tu = t n = t. It has two complex stationary solutions 



± -t± VW^l. (8) 



The equation in variations for (JJJ) is 



Therefore, for < t < 1 both solutions (jHJ are asymptotically stable, while for t > 1 
the solution /i + is asymptotically stable, and fj,~ unstable. A solution of the equation 
with frozen t v approaches its limit exponentially fast, so the value of v near n doesn't 
change significantly while the stabilization occurs. Since by (j2J 

a(n; k) = k (-fj, n + t), 

and the attractor fi + is unique in the case t > 1, we immediately obtain © in that 
case. In the case < t < 1, the solution \i v approaches yT with negative imaginary 
part, because the initial value ~ — h cf. ©• 

Note. The Lemma holds for t = 1 due to the asymptotics derived from [H] (9.3.31-34) 
dkH^\k) 6l/3 (l + iV3)r(2/3) fe ^ ^ 



^ 1} (fc) (i-n/3)r(i/3) 
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Disc of arbitrary radius 



Let u(r, (ft; R, k) be a solution of Problem (1) with wavenumber k outside a circle of 
radius R. Then u{r/R,(ft; l,kR) is a solution of Problem (1) with wavenumber kR 
outside the unit circle. The Dirichlet data for the two functions (as functions of (ft) are 
identical, fn,k{<ft) = (</>)• The Neumann data are related via 

9R,k{<t>) = d r u(r,(ft; R,k)\ r=R = R' 1 gx^R^)- 

Correspondingly, the symbol of the operator M for the disk of radius R is 

a R ( n , k) = R^axin, kR), (9) 

so the limit formula © of Lemma holds with n = n(k,t) = [kRt]. Equivalently, we 
can write the argument of the limit function <7n m (t) as 

n 2tt n , . 

f = M = TP < 10 » 

where L = 2irR is the circumference of the boundary. Notice that the factor 2ir/L is 
the Jacobian d(ft/ds of the boundary parameter change from the arclength s to (ft. 

In the limit R — > oo the disk becomes a half-plane and an analog of the asymptotic 
formula © is an exact formula (|T2|) below. 



Half-plane 



For the Helmholtz equation in the half-plane {x £ R, y > 0), Sommerfeld's radiation 
condition is replaced by a condition that explicitly specifies allowed harmonics in the 
decomposition of any outgoing solution. Namely, two differently behaved families of 
elementary outgoing solutions are given by 

!exp{ia;£ + iy^Jk 2 - £ 2 }, -1 < f < 1, 
exp(^£) exp(-yV£ 2 = k 2 ), 1^1 > L 
The general outgoing solution has the form 

/oo 
f(£)w(x,y;0dt, (11) 

(we don't discuss possible classes to which the function /(£) may belong). 

It is readily seen that /(£) is the Fourier transform of the Dirichlet boundary data 
f(x) = u(x, 0). Differentiating (|11|) with respect to y, we obtain the Fourier representa- 
tion for the Neumann data g(x). The formula for the Dirichlet-to-Neumann operator, 
an analog of reads 

Nf(x) = — / a(0m^d^ 

where the symbol <7(£) = cr(t;;k) is d y w(x,y; £)/w(x,y; ^) 1^=0, i-e. 

a{£;k) = kax^ti/k). (12) 
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Periodic pseudodifferential operators 

The reader can probably see what conclusion we are about to draw from the above 
examples. Let us complete technical preparations, then formulate the main conjecture. 

Recall briefly and informally some basic notions regarding pseudodifferential oper- 
ators on the unit circle S. See JJ| f° r a f un account of the topic. 

Let a((f>, n) be a function on S x Z, which satisfies certain regularity conditions. 
The function a((f>, n) is the discrete symbol of the periodic pseudodifferential operator 
(PPDO) A defined by the formula 

Af(</>) = 53o(^n)/(n)e^. 

neZ 

Here /(</>) is a 27r-periodic function and f(n) its Fourier coefficients. 

The symbol <r(n) introduced in (J3J) does not depend on <j>. Such symbols are called 
constant symbols, and corresponding operators are shift invariant. 

The theory of PPDO applies not only to operators on the unit circle, but to oper- 
ators on any smooth closed curve, since functions on closed curved can be identified 
with 2-7r-periodic functions by reparametrization. 

The symbol a(4>, n) of a PPDO A typically has an asymptotic expansion in decreas- 
ing powers of n. The principal symbol is the leading term in the asymptotics 

a(cj>, ±|n|) = a (0) \n\ a + o(|n| Q ), \n\ — ► oo, 

and a is the order of A. For example, for any domain the operator M is a PPDO of 
order of 1, and for the unit disk its principal symbol is — \n\, cf. 

Theory of PPDO is somewhat simpler than the general theory of pseudodifferential 
operators on compact manifolds (see e.g. The definition of a general PDO uses 

partition of unity. Only the principal symbol can be defined globally. 

The discrete symbol a(4>, n) of a classical PPDO agrees on Z with a symbol a(</>, £) 
defined in the general theory, modulo a function with asymptotics 0(\n\~°°). 

Reconstruction of an operator by its symbol in the general theory assumes that 
operators with smooth kernels are neglected. It isn't convenient when one studies 
double asymptotics (in £ and k), since the behaviour of the neglected part with respect 
to k is not controlled. Correspondence between operators and symbols in the theory of 
PPDO with discrete symbols is strict and preserves full information in both directions. 

Limit Shape Conjecture 

We return to Problem 1 with general boundary T. Denote the length of T by L. Let s 
be the arclength parameter on T (with an arbitrarily chosen starting point), and set 

27T 

ip = s — , <ip <2ir. 
Li 

Consider the operator M as a PPDO (with respect to the parametrization by ip). 
Denote its symbol as ar(i/j,n; k), emphasizing dependence on the wavenumber k. 
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Conjecture. For any fixed t € R and n = n(fc, t) 

o- r (ip,n; k) 



lim 

k— >oo 



[(L/27r)/ct], i/iere exists 

Clim(i)- 



uniformly w.r.t. tp. The universal function o~u m (t) is defined in 
Let us say less formally: 



(Tr(ip,n; k) A; cth 



2tt n 



We can make the conjecture even more readable at the expense of precise mathematical 
meaning. Let us ignore problems associated with definition of a global symbol of 
PDO in the standard theory, where the frequency argument is continuous. Assume 
that <7r(s,£; k) is the symbol of the operator J\f corresponding to the arclength 
parametrization of the boundary. Then 



<r r (s,£; k) 



i^jw^e, 



i<k 
i>k. 



(13) 



Thus the symbol for any boundary parametrized by the arclength is asymptotically 
equal to the exact symbol for the half-plane. This conclusion is hardly surprising given 
that at high frequencies the diffraction process is well localized and (|13|) takes place for 
any disc — see 0, (|10|) — and doesn't refer to curvature. 

Our conjecture has no problems with tangent rays and shadow regions since the 
formula doesn't depend on the boundary data. In particular — in the case of a plane 
incident wave — the direction of incidence has no effect on our claim. One can argue 
that the conjecture has no backing in the case of non-convex scatterers. In that case it 
is supported by numerical results; see the last section of the paper. 

Kirchhoff's approximation 

A relation between the boundary data 
/ and g of an outgoing solution can be 
described alternatively by the impedance 
function n = g/f. It depends on the 
solution. However, according to Kirch- 
hoff's approximation, at high frequen- 
cies the impedance function approaches 
an universal function that depends only 
on the boundary shape. Let us "derive" 
this approximation from the Conjecture. 
Consider an incident plane wave u mc with the wave vector fcko, ||ko|| = 1. Let n 
be the unit normal vector to the boundary V at the given point P E T. Denote by 
9 the angle between ko and n (Fig. 1). The incident wave length A = 2n/k is the 
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distance between wave fronts with equal phases. The boundary value u\nc lr oscillates 
with period A = A/ sin# near the point F. We say that local frequency of iti nc |r at F 
is £ = 2ir/A = ksinO. Assuming Dirichlet's condition for the total field u; nc + u the 
boundary value f = ur also oscillates with local frequency £ at P. 

From a physical point of view, the action of the operator Af amounts to multi- 
plication of local Fourier harmonics by the values of the symbol or at corresponding 
space-frequency locations. In the present case, where the harmonic with frequency £ 
dominates at point P, formula (|13|) implies 

Afu(P) » <7 r (F,£; k)u(P) « iyJk 2 -£ 2 u(P) = ik cos9u(P). (14) 

Fig. 1 shows an illuminated region of the boundary, but the argument holds for a 
shadow region as well. Formula (|14jl can be written in the form 

n(P)nik |(k ,n(P))|, 

which is the classical Kirchhoff approximation formula [§]. A rigorous mathematical 
treatment of Kirchhoff 's approximation (for convex domains) is given in j!2l Ch. X]. 

Insufficiency of the naive local frequency analysis 

The simplistic understanding of the sym- 
bol via local frequencies fails in the fol- 
lowing example. Consider the horse- 
shoe domain $7 as shown on Fig. 2. Let 
two solutions and of Problem 
(|T|) be defined outside £1 as cylindrical 
waves generated by the fictitious sources 
at the points Sj, j = 1,2, inside 
From asymptotics of Hankel's function 
H^\kr) we see that if k\SjP\ > 1, 
then the two solutions yield opposite 
impedances 771(F) ps —772(F) ps — ik. 
The local tangential frequency at F is Fig. 2 

close to for both solutions. Thus it is 

impossible to determine the value or(F, 0; k) consistently by this approach. 
Numerical verification of the Conjecture 

The following algorithm has been used to retrieve the symbol of the operator Af. 
Assume k is given. The algorithm has three free parameters: number of nodes N 
(taken in the form N = 2 m for convenience), and coordinates (xs, Us) of a fictitious 
source inside the domain f2. 
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Algorithm. 

1. Find an equidistant partition of T by N nodes Pj. 

2. Boundary data will be taken from the sample outgoing solution 

u(p) = H^\k\ps\), p$n, 

where S = (xs, Us) is the "source", and P is an observation point. Compute the 
boundary data U = u(Pi), g t = d n u(Pi), i = 1, . . . , N. 

3. Compute discrete Fourier transforms f(n), g(n), n = 0, ... ,N — 1, of the arrays 
{ft}, {gi} using FFT algorithm. Only the first n max Fourier coefficients are considered 
reliable and are used in the sequel. 

4. Find the truncated symbol of a shift-invariant operator that takes / to g: 

cr(n) = g(n)/f(n), n = 0, . . . , n max - 1. 

5. To verify the Conjecture, compare the values k -1 a(n) to an m (2Tvn/kL), where L 
is the length of T. 

We present results obtained for the kite domain p. 70] shown on Fig. 3 and defined 
by the parametric equations 

x(t) = cost + 0.65 cos 2t - 0.65, y(t) = 1.5 sini, t = 0...2vr. 
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1 -0.5 
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Fig. 3: Test domain ("kite") 

The parameters are: k = 200, N = 2 20 , S(— .7, .5). The width of the triangle on Fig. 3 
is equal to 10 wavelengths. In this example, length L = 9.32402 and kL/2Tr ~ 297. 

On Fig. 4, the horizontal coordinate is t = 2irn/kL. Thick lines show the normalized 
real (a), with negative sign, and imaginary (b) parts of the computed approximate 
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symbol, k a(n). Thin lines are the conjectured limit shapes. The true symbol or 
in this case is non-constant, so the approximation by a shift-invariant symbol depends 
on the chosen position of the source. For a source closer to the center of the kite, 
oscillations near t = 1 become smaller. However, in that case the computed values 
near t = 2 oscillate wildly, because corresponding Fourier coefficients f(n) become 
evanescent. 




2.0 t 



Fig. 4: Computed symbol k (?{n) vs <J\i m (t), t = f^n 

The upper bound i max 2.3 on the graphs corresponds to n max = 700 set in the 
computer program. Stabilization of the Fourier coefficients at the upper end of this 
range occurs for the order of discretization N > 2 18 . Obtaining stable values of the 
approximate symbol at larger values of t requires use of larger values of N that grow, 
roughly, exponentially with t. 

A program used for these calculations had a 12 byte long type for floating point 
operations (long double in C). The results obtained with a 8 byte long arithmetics 
(C's type double) were nearly identical. So in the considered example numerical errors 
due to a limited precision are not an issue. 

Conclusion 

The main result is the proposed Limit Shape Formula (|13|) for the symbol of the 
Dirichlet-to-Neumann operator for the standard 2D diffraction problem Q with smooth 
boundary. This asymptotics is independent of the boundary data, of the boundary cur- 
vature, and of convexity assumptions. The limit function cr\{ m {t) defined in (jUJ) varies 
slowly in its argument t ~ const n/k, except near t = 1. These features make the ap- 
proximation (|13fl useful for numerical completion of the boundary data set (u\r, d n u\r), 
which yields the solution u and the radiation pattern by Green's formula. This ap- 
proach includes and supersedes the classical Kirchhoff approximation. We believe that 
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the asymptotics can be enhanced and next, curvature-dependent, term(s) can be found 
from the theory of pseudodifferential operators. In the especially important region, a 
narrow neighbourhood of t = 1, methods for a field near a caustic [I] can be used. 
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